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I.  Motivation  and  Overall  Goals 

The  current  interest  in  the  modeling  and  design  of  emerging  technologies  such  as  very  low 
observable  vehicles,  ground/foliage  penetrating  radars  and  phase  sensitive  components,  imposes 
requirements  on  the  accuracy  and  performance  of  the  computational  tools  far  beyond  the  capa¬ 
bilities  of  existing  techniques.  Furthermore,  the  need  to  accurately  model  the  interaction  of  very 
broad  band  signals  with  electrically  large  and  geometrically  complex  objects,  often  including 
regions  of  advanced  materials,  suggests  that  new  approaches  beyond  integral  equation  based 
approaches  to  electromagnetic  modeling  and  design  be  sought. 

As  has  been  realized  over  the  last  few  years,  the  widely  used  finite-difference  time-domain 
(FD-TD)  method  for  the  time-domain  solution  of  Maxwells  equations  has  a  number  of  very 
unfortunate  properties.  On  one  hand,  the  theoretical  second  order  spatial  and  temporal  accu¬ 
racy  requires  the  use  of  a  significant  number  of  grid  points  per  characteristic  wave  length  to 
accurately  resolve  the  wave  propagation  and  reflection  and  refraction  at  interfaces.  This  again 
translates  into  a  very  large  number  of  degrees  of  freedom  and  makes  the  modeling  of  large  scale 
scattering  and  penetration  problems  prohibitive.  Furthermore,  the  straightforward  use  of  a  sim¬ 
ple  Cartesian  grid  structure  implies  that  material  interfaces  and  metallic  boundaries  be  forced 
to  align  with  the  underlying  grid  structure,  hence  prohibiting  the  accurate  representation  of 
curved  interfaces/boundaries  and  reducing  the  accuracy  of  the  scheme  to  first  order.  A  final, 
and  often  overlooked,  problem  is  the  inability  of  the  FD-TD  approach  to  enforce  the  correct 
jump  conditions  of  the  electric  and  magnetic  field  components  across  material  interfaces.  A ? 
can  be  shown  through  simple  analysis,  this  reduces  the  accuracy  of  the  scheme  to  less  that  first 
order  and  may  even  result  in  nonconvergent  behavior. 

It  is  to  resolve  these  critical,  and  crippling,  issues  that  we  during  the  last  year  have  initiated 
the  development  of  a  new  generation  of  high-order/spectral  schemes  for  the  time-domain  solution 
of  Maxwells  equations  employing  a  truly  unstructured  grid  volume  representation. 

It  is  well  known  that  the  key  technique  to  overcome  the  numerical  errors  associated  the 
classical  FT-TD  approach  is  to  use  high-order/spectral  accuracy  methods,  which  allows  for  a 
dramatic  reduction  of  the  number  of  grid  points  per  characteristic  wavelength  without  sacrificing 
the  overall  accuracy,  i.e.,  for  a  high-order  scheme  one  can  typically  accurately  model  the  wave 
phenomena  with  4-6  points  per  wavelength.  Geometric  flexibility,  hence  overcoming  the  problems 
with  representing  curvilinear  bodies  and  interfaces,  is  ensured  by  using  a  fully  unstructured 
body-conforming  grid  which  also  allows  for  enforcing  the  physically  correct  jump  conditions  on 
the  individual  field  components  across  material  interfaces  and  at  metallic  boundaries. 

The  unstructured  grid  scheme  is  entirely  local  and  information  is  only  exchanged  between 
faces  of  elements,  i.e.,  the  formulation  allows  for  very  efficient  parallel  implementations,  essential 
to  enable  the  modeling  of  large  complex  objects.  Furthermore,  the  use  of  an  standard  unstruc¬ 
tured  grid  allows  for  the  tight  integration  with  existing  key  elements  of  a  design  and  analysis 
loop  such  as  grid  generation  technology,  CAD  systems  and  pre-  and  post-processing  analysis 
tools.  This  latter  point  is  often  underestimated  but  plays  a  key  role  if  a  transition  to  a  larger . 
user  base  is  to  be  successful  and  has  therefore  been  a  guiding  light  throughout  the  development, 
initial  implementation  and  verification. 

II.  Current  State  of  Development  and  Main  Achievements 

In  the  following  we  shall  review  the  main  accomplishments  during  the  period  of  the  project 
as  related  to  the  continued  development  of  high-order/spectral  time-domain  methods.  We  shall 
overview  the  status  of  the  development  of  structured  grid  spectral  methods,  the  more  recent 
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Figure  1:  Left:  Typical  multi-domain  grid  for  the  solution  of  the  scattering  by  a  missile.  Right:  RCS(#,0) 
for  a  missile  subject  to  axial  illumination  by  a  horizontally  polarized  plane  wave.  The  reference  solution, 
marked  by  was  obtained  by  a  method-of-moments  technique. 

development  of  high-order/spectral  methods  on  unstructured  grids,  the  main  results  of  the 
continued  efforts  in  the  analysis  of  PML  methods,  as  well  as  the  adaptation  of  the  boundary 
variation  methods  to  problems  in  diffractive  optics.  We  shall  conclude  by  discussing  efforts 
related  to  finite-difference  time-domain  methods  and  techniques  to  improve  on  the  widely  used 
Yee  scheme  to  enable  correct  treatment  of  interfaces  and  boundaries. 

II.  1  Structured  Grid  Multi-Domain  High-Order/Spectral  Methods 

The  straightforward  application  of  spectral  methods  in  several  space  dimensions  requires 
that  the  computational  domain  is  a  square  (in  two  dimensions)  or  a  cube  (in  three  dimensions) 
with  the  grid  predefined  by  a  tensorproduct  of  Gaussian  points  to  assure  the  spectral  accuracy. 

In  the  structured  grid  multi-domain  formulation  these  squares  and  cubes  take  the  role  of 
the  fundamental  building  blocks  by  which  the  computational  domain  is  filled  to  enable  a  high- 
order /spectral  solution  of  Maxwells  equations  in  a  geometrically  complex  domain.  Filling  the 
computational  domain  is  done  in  a  completely  body-conforming  way,  representing  curvilinear 
material  interfaces  and  metallic  boundaries  with  an  accuracy  equal  to  that  of  the  approximation. 
The  computations  in  these  elements  involve  the  approximations  of  equations  obtained  by  means 
of  curvilinear  transformations  of  Maxwells  equations. 

To  connect  the  domains  we  exploit  the  hyperbolic  nature  of  Maxwell  equations  along  with 
the  uniquely  defined  outward  pointing  normal  vectors  at  the  surfaces  of  the  curvilinear  elements 
to  uniquely  determine  which  information  leaves  and  which  enters  the  domain.  Continuity  of 
the  characteristic  variables  is  enforced  strongly  in  a  way  entirely  consistent  with  the  physics 
of  the  problem.  For  interfaces  placed  between  two  materials  we  enforce  the  electromagnetic 
jump  conditions  to  ensure  correct  behavior  of  the  fields  along  such  interfaces.  This  approach 
is  likewise  employed  at  geometric  singularities  such  as  vertices  and  edges  combined  with  very 
weak  filtering  to  ensure  stability  of  the  approximation. 

Using  these  ideas,  we  have  demonstrated  the  prospects  of  using  high-order/spectral  multi¬ 
domain  methods  for  the  accurate  and  efficient  solution  of  nontrivial  problems  in  electromagnetics 
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Figure  2:  Numerical  RCS  result  of  scattering  by  a  sphere  (left)  and  a  finite  cylinder  (right)  of  a  900 MHz 
plane  wave.  Both  bodies  are  composed  of  lossy  di-electric  materials  with  er  =  43.0,  pr  =  1.0,  and 
a  =  0.8Sm-1.  The  sphere  has  a  radius  of  0.1  meter,  i.e.,  ka  =  0.67T,  and  the  length  of  the  cylinder  is 
0.67m.  The  dashed  line  corresponds  to  the  Mie-series  result  (left)  and  a  finite  element  solution  (right) 
respectively. 


involving  scattering  by,  and  penetration  into,  geometrically  complex  two-  and  three-dimensional 
bodies. 

As  an  example  of  current  capabilities,  consider  the  problem  of  scattering  from  an  axisymmet- 
ric  missile,  assumed  to  be  perfectly  conducting.  In  Fig.  1  we  illustrate  a  typical  structured  grid 
for  the  computation,  emphasizing  the  fully  body  conforming  grid  and  the  geometric  flexibility 
of  the  method.  The  validity  of  the  results  are  confirmed  by  a  direct  comparison  with  a  standard 
method-of-moments  computations. 

For  the  modeling  of  problems  containing  advanced  materials,  the  development  is  currently 
ongoing  and  much  additional  work  is  needed.  For  problems  with  piecewise  homogeneous  ma¬ 
terials,  both  two  and  three-dimensional  tests  have  been  successfully  completed,  confirming  the 
accuracy  and  overall  performance  for  scattering  by  spheres,  cylinders,  cubes  etc. 

As  an  example  of  the  modeling  of  scattering  by  lossy  materials  we  show  in  Fig.  2  the  RCS 
computed  using  the  three-dimensional  framework  for  monochromatic  scattering  by  a  very  dense 
sphere  and  a  very  dense  cylinder  with  the  latter  being  compared  with  results  obtained  using 
a  very  highly  resolved  finite  element  approach.  The  parameters  of  the  materials  are  chosen  to 
assemble  those  of  human  tissue. 

Both  the  two  and  the  three  dimensional  codes  exploit  the  use  of  Perfectly  Matched  Layers 
(PML)  to  terminate  the  computational  domains  and  specially  designed  PML  methods  for  high- 
order  methods  have  been  implemented  and  thoroughly  tested. 

To  enable  the  computational  modeling  of  electrically  large  problems  involving  advanced 
materials,  the  potential  for  efficient  implementations  is  critical.  The  existing  codes  have  all  been 
implemented  using  MPI  to  support  fully  parallel,  distributed  memory  execution  at  a  number  of 
platforms,  showing  superior  parallel  performance 

These  results  clearly  establish  that  the  recently  developed  two-  and  three-dimensional  struc¬ 
tured  grid  high-order/spectral  multi-domain  schemes  are  geometrically  flexible,  capable  of  han¬ 
dling  real  materials  and  allow  for  efficient  implementations  to  enable  the  modeling  of  large  and 
complex  realistic  configurations. 
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II. 2  Unstructured  Grid  Multi-Domain  High-Order/Spectral  Methods 

i A.  concentrated  effort  has  recently  been  initiated  to  develop  methods  for  solving  the  time 
domain  Maxwells  equations  in  complex  geometries  and  media  using  an  fully  unstructured  grid. 
The  USEMe  code  (Unstructured  Sprectral  Element  Method)  is  the  resulting  implementation 
of  a  numerical  scheme  which  employs  nodal  high-order /spectral  multi-domain  methods.  Each 
component  of  the  code  has  been  carefully  designed  with  efficiency,  accuracy,  flexibility,  and 
robustness  in  mind. 

The  standard  geometric  building  block  by  which  the  computational  domain  is  gridded  in 
USEMe  is  the  triangular  element  in  2D,  and  the  tetrahedral  element  in  3D.  This  allows  for  the 
use  of  standard  finite  element  or  finite  volume  type  meshes  using  these  elements  to  describe 
arbitrarily  complex  domains,  hence  overcoming  one  of  the  main  concerns  associated  with  the 
structured  grid  methods.  The  use  of  triangles  and  tetrahedrons  remains  of  the  most  important 
design  principles  of  the  method  as  it  enables  a  tight  integration  with  industry  standard  mesh 
generators. 

USEMe  is  based  on  penalty/discontinuous  Galerkin(DGM)  type  scheme,  which  is  a  particular 
form  of  a  penalty  scheme  that  guarantees  elemental  conservation  and  global  stability.  Neighbor¬ 
ing  elements  are  patched  together  through  the  use  of  penalty  terms  to  account  for  jumps  across 
their  shared  boundary.  ThU  formulation  has  the  advantage  that  it  allows  material  coefficients 
to  vary  discontinuously  at  the  interface  between  two  elements,  while  still  preserving  high  order 
accuracy.  Additionally  chis  allows  for  efficient  communication  in  parallel  computations. 

A  thorough  convergence  analysis  for  the  proposed  framework  has  been  completed,  putting  the 
expectations  on  firm  ground  and  verifying  that  one  can  indeed  expect  the  properties  speculated 
during  the  initial  developments  and  essential  as  a  motivation  for  the  chosen  approach.  The 
analysis  also  includes  bounds  on  the  divergence  error,  confirming  that  the  computed  results 
remain  divergence  free  to  the  order  of  the  scheme  away  from  singular  points. 

A  general  two-dimensional  unstructured  grid  code,  based  on  the  use  of  curvilinear  triangles, 
has  been  implemented  and  tested  extensively,  including  by  third  party  users.  The  code  is  flexible, 
versatile,  and  robust,  allowing  for  the  modeling  a  very  general  two-dimensional  (TE  and  TM) 
time-domain  scattering,  penetration  and  radiation  problems  involving  di-electric  and  magnetic 
materials,  possibly  of  a  lossy  nature. 

The  discretization  provides  a  fully  body-conforming  high  order  representation  of  a  mate¬ 
rial  scatterer,  hence  overcoming  the  staircasing  phenomena  evident  in  common  finite  difference 
schemes.  In  Figure  3  we  illustrate  some  elements  which  have  an  edge  on  the  boundary.  Their 
nodes  have  been  mapped  to  the  boundary  and  the  internal  nodes  blended  linearly  from  the 
boundary.  The  full  triangulated  region  is  likewise  illustrated. 

As  an  example  of  the  expected  accuracy,  we  compute  plane  wave  scattering  from  a  PEC 
cylinder  with  the  non-uniform  near- field  mesh  illustrated  in  Fig.  3.  This  is  an  excellent  test 
case  as  it  allows  ns  to  determine  the  performance  of  the  PEC  boundary  conditions,  PML  bound¬ 
ary  truncation,  RCS  calculation,  curvilinear  body  representation  and  high  order  approximation 
techniques.  Comparing  the  computed  z  component  of  the  electric  field  with  the  exact  solution 
over  the  interior  domain,  Table  1  confirms  exponential  convergence.  Numerous  other  results  for 
scattering  and  penetration  into  simple  simple  homogeneous  scatters  have  shown  similar  accu¬ 
racy.  Other  more  complex  problems  involving  materials  and  singular  solution  behavior  display 
a  very  similar  accuracy  and  overall  performance. 

Serial  as  well  as  fully  parallel  implementations  has  been  completed  and  various  tools  such 
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Figure  3:  Left:  Seventh  order  boundary  fitted  triangular  elements.  Right:  Full  triangulated  region  for 
two-dimensional  scattering  with  PML. 


Expansion 

Loo  error 

8 

3.2E-02 

10 

2.1E-03 

12 

6.2E-04 

Table  1:  L^o  error  over  the  interior  domain  as  function  of  expansion  order  for  scattering  of  a  TM  plane 
wave  from  a  PEC  cylinder  (ka  =  77r). 


as  near-to-far-field  transformations  and  highly  efficient  absorbing  layer  techniques  have  been 
implemented  and  tested.  The  code  has  been  interfaced  with  a  public  domain  grid-generator  to 
enable  simple  and  user-friendly  grid-generation. 

The  development  of  a  general  three-dimensional  framework  is,  by  the  sheer  complexity  of  the 
problem,  a  daunting  task  and  it  remains  a  work  in  progress.  A  general,  three-dimensional  version 
of  the  scheme  has  been  implemented  and  initial  tests  has  been  conducted,  verifying  the  expected 
accuracy  and  overall  performance.  So  far,  these  test  cases  has  been  limited  to  problems  involving 
purely  metallic  scatters,  although  this  remains  less  of  a  concern  as  pure  metallic  scattering  often 
provides  the  most  challenging  test  cases. 

As  a  simple  test  of  the  current  capability  of  the  framework  we  show  in  Fig.  4  the  bistatic 
RCS  for  a  ka  =  10  PEC  sphere  as  computed  using  the  unstructured  framework  and  compared 
with  the  analytic  Mie-series.  As  expected  we  find  excellent  agreement. 

As  a  more  eye-catching,  and  preliminary,  result  Figure  5  demonstrates  that  the  code  is 
currently  capable  of  simulating  electromagnetic  plane  wave  scattering  from  a  complex  F15  type 
geometry.  Even  though  this  is  a  preliminary  result,  the  code  is  capable  of  representing  much  more 
complex  media  and  also  more  complex  and  higher  frequency  signals  than  this  monochromatic 
example.  This  result  was  obtained  by  running  on  a  single  processor  workstation  using  third 
order  accurate  elements.  Higher  resolution  calculations  are  being  performed  in  parallel  on  an 
IBM  SP. 
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Figure  4:  Bistatic  RCS  for  ka  =  10  as  computed  using  the  unstructured  grid  framework  with  an  6’th 
order  basis. 

A  fully  parallel  version,  implemented  in  MPI  and  Fortran/C,  has  been  tested,  showing  ex¬ 
cellent  parallel  speed  up  which  remains  a  crucial  property  to  enable  the  modeling  of  large  scale 
problems  of  realistic  complexity.  The  parallel  performance  is  illustrated  in  Table  2,  confirming 
that  the  code  attains  up  to  95  percent  parallel  efficiency  as  the  number  of  processors  is  increased. 


Expansion 

Total  Degrees 

#  of  processors 

Order 

of  Freedom 

2 

4 

8 

16 

3 

7,380,000 

5.6 

3.4 

1.8 

.95 

4 

14,760,000 

** 

7.4 

3.9 

2.3 

5 

25,830,000 

** 

** 

7.8 

4.3 

Table  2:  Wall  clock  time  per  Runge  Kutta  integration  stage.  Timings  taken  on  an  IBM  SP2  at  CASCV, 
Brown  University.  (**  implies  that  there  was  not  enough  memory  local  to  the  nodes) 

The  code  has  been  interfaced  with  Gambit/Fluent  grid-generator  to  allow  for  a  user-friendly 
approach  to  the  daunting  of  three-dimensional  grid-generation.  This  interfacing  has  been  accom¬ 
plished  through  a  pre-processing  module  that  enables  the  use  of  a  standard  finite-element /finite- 
volume  grid  in  a  fully  body-conforming  higli-order  accuracy  framework.  It  is  worth  emphasizing 
that  there  is  nothing  special  about  the  chosen  grid-generator  and  any  finite-element /finite- 
volume  grid  generator  can  be  used  as  a  frontend  to  USEMe,  provided  only  that  a  simply  format- 
interface  be  implemented. 

II. 3  Absorbing  Boundary  Conditions  for  Time-Domain  CEM 

One  of  the  central  issues  in  time  domain  CEM  is  the  development  of  techniques  for  the 
truncation  of  the  infinite  computational  domain  in  such  a  way  that  outgoing  waves  leave  the  do¬ 
main  without  reflections  which  could  otherwise  reenter  and  eventually  falsify  the  computational 
results. 
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Figure  5:  Nose  to  tail  component  of  the  reflected  magnetic  field  from  an  incident  plane  wave. 


In  1994  Berenger  introduced  a  new  methodology  for  designing  such  schemes  by  the  introduc¬ 
tion  of  an  absorbing  layer,  to  which  he  attributed  material  properties  that  modify  the  vacuum- 
equations  so  that  the  field  strength  decays.  In  his  original  paper  these  material  properties  were 
based  on  a  mathematical  construct,  involving  splitting  the  transverse  magnetic  field. 

We  have  previously  conducted  a  detailed  mathematical  analysis  of  the  PML  formulation 
applied  to  the  two  dimensional  transverse-electric  mode  (TE)  of  Maxwell’s  equations.  The 
conclusions  of  this  analysis  hold  also  for  all  other  forms  of  the  equations.  We  were  interested 
in  the  well-posedness  of  the  PML  formulation,  and  therefore  investigated  the  pure-initial  value 
problem.  The  main  conclusion  is  that  the  PML  split-form  of  Maxwell’s  equations  is  only  weakly 
well-posed  and  therefore  its  solution  could  diverge  under  small  perturbations.  This  theoretical 
result  has  some  far  reaching  conclusions.  First,  it  indicates  the  the  split  field  PML  will  have 
problems  as  the  required  number  of  points  in  the  PML  layer  will  increase  with  time.  Moreover, 
since  all  physical  system  describing  reality  must  be  strongly  well  posed,  the  PML  can  not 
represent  any  physical  layer. 

In  the  case  of  the  transverse  electric  (TE)  case  or  transverse  magnetic  (TM)  these  consider¬ 
ations  lead  to  a  set  of  4  p.d.e.’s,  different  from  those  given  by  Berenger. 

The  dimensionless  form  of  the  equations  based  on  the  Lorentz- material  model,  is: 

0EX  _  8H  8Ey  _  _0H  _ 

dt  dy  '  dt  dx  y 
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OH__dE1^dEy_  dJ  _  dH 
dt  dy  dx  °  ’  dt  °  dy 

where  J  is  a  polarization  current. 

We  showed  that  the  set  of  4  p.d.e.’s  can  be  reduced  to  a  set  of  3  inhomogeneous  Maxwell 
equations  augmented  by  a  temporal  O.D.E.  This  means  that  the  original  set  of  equations  in  the 
layer  is  strongly  well  posed.  We  then  provide  a  closed  form  solution  for  plane  -waves  traveling 
in  the  absorbing  layer,  under  the  assumption  of  periodicity  in  y,  for  2  cases  -  the  semi-infinite 
layer  and  the  finite  layer.  The  analysis  confirmed  that  the  infinite  layer  is  a  true  PML  but  also 
exposed  that  the  finite  width  Lorentz  material  PML  may  exhibit  locally  growing  solutions  before 
they  finally  decay. 

The  method  of  deriving  absorbing  layers  from  physical  considerations  has  limitations.  In 
particular  the  absorbing  properties  of  the  layers  can  not  be  controlled.  We  have  developed  a 
mathematical  methodology  for  the  construction  of  PMLs.  This  methodology  is  more  flexible 
than  the  physically  motivated  one  and  it  allows  control  of  the  decay  properties  of  the  solution 
in  the  absorbing  layer.  One  of  the  possible  sets  is  given  by 
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The  solution  here  is  bounded  by  a  function  that  does  not  grow  in  any  direction 


The  success  of  the  PML  method  is  crucial  in  applying  high  order  accuracy  methods  to  the 
simulations  of  Maxwell  equations.  After  all,  there  is  no  point  in  achieving  high  accuracy  in  the 
computational  domain  if  the  reflections  introduce  large  errors.  The  PML  methodology  minimize 
the  numerical  reflections. 


However  it  was  found  that  all  the  known  PML  allow  linear  (and  non-physical)  growth  in  time 
when  simulating  pulses.  It  seems  that  the  growth  ig  triggered  when  the  pulse  leaves  the  domain 
and  the  solution  vector  approach  constant  in  space.  A  simple  analysis  explains  the  situation. 
Consider  the  physically  based  PML,  and  assume  that  the  spatial  derivative  vanish,  then  the 
equations  for  Ex  and  J  decouple  and  they  display  a  linear  growth  in  time.  It  is  easy,  based  on 
this  analysis,  to  find  a  ’’cure”  for  this  problem.  However  it  is  not  clear  if  the  corrected  method  is 
still  a  PML.  This  problem  remains  open.  The  numerical  experiments  indicate  that  the  corrected 
.physically  motivated.  PML  maintains  its  absorption  properties  even  with  the  time  stabilizing 
term.  A  better  analysis  and  more  experiments  are  needed  to  have  the  full  understanding  of  this 
very  important  question. 
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II. 4  Boundary  Variation  Techniques  for  Waveguide  Holograms 

It  has  recently  been  established  by  Bruno  and  Reitich  that  solutions  to  electromagnetic 
diffraction  by  a  periodic  structure  depend  analytically  on  the  variations  of  the  interface.  In  other 
words,  diffraction  from  a  periodic  grating  can  be  determined  from  knowledge  of  reflection  and 
refraction  at  a  plane  interface  and  the  diffraction  problem  be  solved  by  analytic  continuation. 
Using  this  result,  a  high-order  perturbation  technique  was  developed  for  finite-size  perturba¬ 
tions  and  successfully  used  it  in  modeling  diffraction  by  two  and  three-dimensional  metallic  and 
transmission  gratings. 

To  utilize  and  further  generalize  on  these  results  we  have  initiated  work  on  the  formulation 
and  implementation  of  boundary  variation  techniques  for  the  analysis  and  modeling  of  waveguide 
grating  couplers.  These  are  diffractive  optical  devices  in  which  a  guided  wave  in  a  thin-film 
waveguide  is  coupled  to  free-space  radiation  through  a  surface  relief.  These  problems  are  vastly 
different  from  what  was  previously  considered.  Among  other  things,  the  illuminating  wave 
is  a  guided  wave  with  evanescent  tails.  Furthermore  the  diffraction  process  includes  multiple 
refraction/reflection  processes.  Finally,  we  have  extended  the  analysis  to  structures  of  finite 
extent. 


M 

P  (TO"2) 

2 

1.208589 

4 

1.220834 

6 

1.220247 

8 

1.220243 

10 

1.220243 

Table  3:  Power  in  the  -1st  diffraction  order  for  different  number  of  terms  [M,M]  in  the  Pade 
approximation  to  the  power  series  expansion. 


As  a  first  example  of  the  performance  of  the  boundary  variation  scheme,  consider  a  waveguide 
structure  consisting  of  a  core  layer  with  refractive  index  n  =  1.45  and  thickness  d\  =  0.8, 
sandwiched  between  two  cladding  layers  of  refractive  index  n—  1.4.  The  top  cladding  layer  has 
a  finite  thickness  of  d2  =  1  and  above  this  layer  is  air  with  n  —  1.  For  the  fundamental  TE 
mode  this  geometry  yields  an  effective  index  of  1.4213. 

We  consider  a  cosine  surface  relief 

f5{z)  =  ,4  cos  (y2) 

Let  us  first  study  the  cbm'srgence  of  the  scheme.  As  we  do  not  have  an  analytic  solution  to 
compare  with,  we  are  unable  to  compute  the  error.  Rather,  we  look  at  the  power  coupled  to 
the  -1st  diffraction  order  as  we  increase  the  number  of  terms  in  the  Pade  approximation  to  the 
Taylor  series  expansion  of  the  global  solution.  Table  3  confirms  the  convergence  in  the  case  of 
A  =  0.1  as  the  number  of  terms  in  the  power  series  expansion  increases. 

Let  us  now  demonstrate  the  use  of  the  proposed  boundary  variation  method  to  study  aperi¬ 
odic  surface-relief  gratings  couplers  of  finite  length.  Clearly,  as  we  use  a  periodic  Rayleigh  series 
expansion  for  the  radiated  fields,  the  grating  is  implicitly  forced  to  be  periodic.  However,  as 
we  shall  demonstrate,  choosing  a  sufficiently  large  total  length  of  the  computational  domain  as 


10 


compared  with  the  length  of  the  finite  surface  relief,  the  results  becomes  consistent  with  those 
obtained  using  a  method  dealing  with  truly  finite  gratings. 

For  the  FGC  surface  relief  we  use  the  generic  profile 

cos[27T  (a0  +  ai  {z  -  2o))  (z  ~  *o)]  (2) 

where  A  is  the  amplitude,  w  is  the  width  of  the  exponentially  truncated  relief,  zq  is  the  center  of 
the  relief,  ao  =  1/A  for  the  unchirped  relief,  and  oi  is  the  chirp  parameter.  We  consider,  as  an 
example,  the  parameters  A  =  0.25,  ao  =  1.4213,  a\  =  0.005,  and  w  =  3  and  show  in  Fig.  6  that 
there  is  an  excellent  agreement  between  the  farfield  patterns  of  the  results  obtained  even  for  this 
relatively  deep  surface  relief  grating  and  those  computed  with  a  rigorous  spectral  multi-domain 
scheme  solving  Maxwells  equations  directions  and  to  very  high  order. 


fs{z)  =  A  exp 


z  -  z0 


w 


Figure  6:  Left:  Farfield  radiation  patterns  for  FGC  computed  using  the  boundary  variation 
method  (dashed)  and  a  spectral  multi-domain  method  (solid).  Right:  Focusing  beam  emanating 
from  grating  coupler  at  x=0. 


AM  f  time 


BV 

BV 

SC 

BV 

SC 

0.1 

7 

47.1 

47.0 

46s 

30h 

'  •  0.2 

13 

46.3 

46.2 

390s 

64h 

0.3 

17 

45.2 

45.1 

1081s 

126h 

Table  4:  Key  figures  for  spectral  collocation  (SC)  and  the  proposed  boundary  variation  (BV) 
computations.  A  is  the  amplitude  of  the  surface  relief.  An  [M,M]  Pade  approximant  is  used 
for  the  BV  method,  f  is  the  focal  length  normalized  with  the  free-space  wavelength,  and  time 
reflects  the  total  computation  time. 
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We  find  excellent  correspondence  in  the  farfield  maintained  in  the  nearfield  even  though 
the  intensity  is  slightly  lower  for  the  BV  method  for  all  amplitudes,  which  may  be  due  to  not 
accounting  for  multiple  reflections.  Looking  at  the  computation  time  in  Table  4,  it  is  evident 
that  the  use  of  the  approximate  boundary  variation  method  certainly  pays  off:  While  we  find 
excellent  agreement  with  the  rigorous  spectral  collocation  method,  we  find  a  reduction  in  the 
computation  time  exceeding  a  factor  of  2000  is  achievable,  a  fact  which  calls  for  the  future  use 
of  the  method  as  the  forward  solver  in  an  optimization  scheme. 

II. 5  Embedded  Interface  Finite-Difference  Methods 

While  it  is  well  known  4th  order  methods  perform  significantly  better  than  2nd  order  methods 
when  solving  wave  problems,  it  is  also  well  known  that  stable  4th  order  finite  difference  schemes 
can  be  constructed  on  equidistant  Cartesian  grids  even  as  one  approaches  the  boundary  and 
one-sided  closures  are  required.  A  central  difficulty,  however,  is  introduced  in  the  need  to  treat 
geometrically  complex  objects.  Rather  than  adapting  a  multi-element  formulation  we  wish  to 
overcome  this  complication  while  maintaining  a  very  simple  Cartesian  grid  and  embed  the  whole 
computational  problem  into  the  grid.  The  success  of  such  embedding  techniques  hinges  critically 
on  the  formulation  of  the  finite  difference  stencils  in  a  way  that  allows  for  the  including  the 
position  of  the  embedded  objects  as  well  as  the  boundary  conditions  of  the  field  on  the  surface 
of  the  objects. 

Attempts  to  construct  stable  and  well  behaved  finite  difference  schemes  that  allows  one  to 
use  a  simple  Cartesian  grid  for  solving  partial  differential  equation  in  geometrically  complex 
settings  is  not  new.  In  particular  for  solving  the  incompressible  Navier-Stokes  equations  can 
one  find  a  wealth  of  methods  which  all  seem  to  have  their  origin  in  the  immersed  boundary 
method.  For  problems  involving  transparent  boundaries,  as  is  generally  the  case  for  problems  in 
electromagnetics,  it  is  necessary  to  seek  an  approach  in  which  one  directly  imposes  the  correct 
interface  conditions  on  the  solution  across  the  interface.  An  approach  used  extensively  in  the 
context  of  electromagnetics  is  to  simply  assign  averaged  material  properties  to  grid  points  at 
the  material  interface,  although  we  have  shown  this  approach  yields  1st  order  accuracy  at  best. 

In  very  recent  work,  we  have  proposed  a  novel  2nd  order  finite  difference  based  approach  for 
solving  Maxwells  equations  but  directly  applicable  to  wave  problems  in  general.  Contrary  to 
previously  proposed  methods,  it  is  equally  applicable  to  the  scalar  and  the  system  case,  the  bulk 
of  the  work  involved  in  treating  the  material  interfaces  is  performed  in  a  preprocessing  phase, 
and  the  scheme  handles  generally  curved  internal  interfaces  as  well  as  fully  reflecting  boundaries 
with  equal  ease  and  in  a  uniform  way. 

Although  the  thorough  theoretical  developments  are  completed  for  the  one-dimensional  prob¬ 
lem  only,  the  approach  generalizes  to  several  dimensions  without  complications.  Consider,  as 
an  example,  the  electromagnetic  scattering  of  a  plane  wave  impinging  on  a  two-dimensional 
di-electric  cylinder.  In  Fig.  7  is  shown  the  temporal  behavior  of  the  global  error  as  well  as  a 
resolution  study  at  a  specific  time.  We  see,  as  expected,  that  incorrect  staircased  treatment  of 
boundaries  and  their  position  severely  limits  the  accuracy  of  the  simple  Yee  scheme.  Contrary 
to  this,  the  new  scheme  is  truly  second  order  and  typically  yields  at  least  an  order  of  magnitude 
improvement  in  accuracy  over  the  Yee  scheme  for  the  same  resolution  and,  thus,  the  same  com¬ 
putational  work.  In  other  words,  one  can  improve  the  accuracy  at  little  or  no  cost  or  one  can 
obtain  results  of  a  comparable  accurate  accuracy  at  dramatically  reduced  computational  cost. 

Current  efforts  include  the  continued  verification  of  the  two-dimensional  framework  for  both 
polarizations  as  well  as  the  initial  phase  of  the  implementation  and  testing  of  a  general  three- 
dimensional  framework. 
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Figure  7:  In  a)  we  show  the  temporal  dependence  of  the  global  L2  error  of  E  for  different  resolu¬ 
tion  in  terms  of  the  free-space  wavelength  for  the  staircased  and  non  staircased  approximation. 
In  b)  is  shown  the  global  error  at  t=25,  illustrating  the  expected  convergence  rate. 

III.  Papers  Citing  Support  from  Grant 
III.l  Appeared  or  Accepted 

•  A.  Ditkowski,  K.  H.  Dridi,  and  J.  S.  Hesthaven,  1999,  Convergent  Cartesian  Grid  Methods 

for  Maxwells  Equations  in  Complex  Geometries  ,  J.  Comput.  Phys.  -  to  appear. 

•  K.  H.  Dridi,  J.  S.  Hesthaven,  and  A.  Ditkowski,  1999,  Staircase  Free  Finite- Difference 

Time-Domain  Formulation  for  General  Materials  in  Complex  Geometries ,  IEEE  Trans. 
Antennas  Propaga.  -  to  appear. 

•  D.  Gottlieb  and  J.  S.  Hesthaven,  2000,  Spectral  Methods  for  Hyperbolic  Problems ,  J.  Comput. 

Appl.  Math.  -  to  appear. 

•  T.  Warburton,  L.  Pavarino  and  J.S.  Hesthaven,  2000  A  Pseudo- Spectral  Scheme  for  the  In¬ 

compressible  Navier-Stokes  Equations  Using  Nodal  Elements,  J.  Comput.  Phys.  163(2000), 

pp.  1-21. 

•  V.  Zharnitsky,  E.  Grenier,  S.  K.  Turitsyn,  C.  K.  R.  T.  Jones,  and  J.  S.  Hesthaven,  2000, 

Ground  States  of  Dispersion  Managed  NLS ,  Phys.  Rev.  E.  62(5),  pp.  7358-7364. 

•  J.  S.  Hesthaven  and  C.  H.  Teng,  2000,  Stable  Spectral  Methods  on  Tetrahedral  Elements , 

SIAM  J.  Sci.  Comput.  21(6),  pp.  2352-2380. 

•  P.  G.  Dinesen  and  J.  S.  Hesthaven,  2000,  A  Fast  and  Accurate  Boundary  Variation  Method 

for  Diffrative  Gratings ,  J.  Opt.  Soc.  Am.  A  17(9),  pp.  1565-1572. 

•  J.  S.  Hesthaven,  2000,  Spectral  Penalty  Methods ,  Appl.  Numer.  Math.  33(1-4),  pp.  23-41. 

•  B.  Yang  and  J.  S.  Hesthaven,  2000,  Multidomain  Pseudo  spectral  Computation  of  MaxweWs 

Equations  in  3-D  General  Curvilinear  Coordinates ,  Appl.  Numer.  Math.  33(1-4),  pp.  281- 
289. 


13 


•  P.  G>  Dinesen,  J.  S.  Hesthaven  and  J.  P.  Lynov,  2000,  A  Pseudo  spectral  Collocation  Time- 

Domain  Method  for  Diffractive  Optics,  Appl.  Numer.  Math.  33(1-4),  pp.  199-206. 

•  J.  S.  Hesthaven,  P.  G.  Dinesen  and  J.  P.  Lynov,  1999,  Spectral  Collocation  Time-Domain 

Modeling  of  Diffractive  Optical  Elements,  J.  Comput.  Phys.  155(1),  pp.  287-306. 

•  S.  Abarbanel,  D.  Gottlieb  and  J.  S.  Hesthaven,  1999,  Wellposed  Perfectly  Matched  Layers 

for  Advective  Acoustics,  J.  Comput.  Phys  154(2),  pp.  266-283. 

•  P.  G.  Dinesen,  J.  S.  Hesthaven,  J.  P.  Lynov  and  L.  Lading,  1999,  Pseudospectral  Method  for 

the  Analysis  of  Diffractive  Optical  Elements,  J.  Opt.  Soc.  Am.  A  16(5),  pp.  1124-1130. 

•  B.  Yang  and  J.  S.  Hesthaven,  1999,  A  Pseudospectral  Method  for  Time-Domain  Computation 

of  Electromagnetic  Scattering  by  Bodies  of  Revolution,  IEEE  Trans.  Antennas  Propaga. 
47(1),  pp.  132-141. 

•  P.  G.  Dinesen  and  J.  S.  Hesthaven,  2000,  Analysis  of  Grating  Couplers  Using  the  Boundary 

Variation  Method.  In  Technical  digest.  Diffractive  optics  and  micro-optics  meeting  and  table 
top  exhibit,  Quebec  City  (CA).  Optical  Society  of  America,  OSA  Technical  Digest  series, 
pp.  84-86. 

•  P.  G.  Dinesen  and  J.  S.  Hesthaven,  2000,  Rigorous  3-D  Analysis  of  Focusing  Grating  Couplers 

Using  a  Spectral  Collocation  Method.  In  Technical  digest.  Diffractive  optics  and  micro¬ 
optics  meeting  and  table  top  exhibit,  Quebec  City  (CA).  Optical  Society  of  America,  OSA 
Technical  Digest  series,  pp.  81-83. 

•  P.  G.  Dinesen,  J.  S.  Hesthaven, and  J.  P.  Lynov,  2000,  Rigorous  Analysis  of  Focusing  Grat¬ 

ing  Couplers  Using  a  Time-Domain  Spectral  Collocation  Method.  In  Diffractive/holographic 
technologies  and  spatial  light  modulators  7.  Optoelectronics  2000,  San  Jose,  CA.  Cindrich,  I.; 
Lee,  S.H.;  Sutherland,  R.L.  (eds.),  (International  Society  for  Optical  Engineering,  Belling¬ 
ham,  WA,  2000),  Proceedings  of  SPIE  3951,  pp.  11-19. 

•  P.  G.  Dinesen,  J.  S.  Hesthaven, and  J.  P.  Lynov,  2000,  Rigorous  Three-Dimensional  Analysis 

of  Surface- Relief  Gratings  Using  a  Spectral  Collocation  Method.  In  Diffractive/holographic 
technologies  and  spatial  light  modulators  7.  Optoelectronics  2000,  San  Jose,  CA.  Cindrich,  I.; 
Lee,  S.H.;  Sutherland,  R.L.  (eds.),  (International  Society  for  Optical  Engineering,  Belling¬ 
ham,  WA,  2000),  Proceedings  of  SPIE  3951,  pp.  2-10. 

•  K.  Dridi  and  J.S.  Hesthaven,  1999,  N-space  Staircase-Free  Finite- Difference  Time-Domain 

Formulation  for  Arbitrary  Material  Distributions:  Numerical  Investigations  on  a  Focusing 
Grating  Coupler  in  Dielectric  Waveguides.  Proc.  of  Integrated  Photonics  Research  IPR  99, 
Santa  Barbara,  CA.  pp.  250-252. 

•  P.  G.  Dinesen,  L.  Lading,  J.  P.  Lynov  and  J.  S.  Hesthaven,  1998,  Waveguides  and  Diffractive 

Elements  for  Non-Contact  Sensors:  Analysis  Proc.  of  Diffractive  Optics  and  Micro-Optics, 
Hawaii,  pp.  209-211. 

•  J.  S.  Hesthaven,  P.  G.  Dinesen  and  J.  P.  Lynov,  1998,  Pseudospectral  Time-Domain  Modeling 

of  Diffractive  Optical  Elements.  Proc.  of  The  14'th  Annual  Review  of  Progress  in  Applied 
Computational  Electromagnetics,  Monterey,  CA.  pp.  858-865. 


14 


•  B.  Yang,  D.  Gottlieb  and  J.  S.  Hesthaven,  1997,  On  the  Use  of  PML  ABC’s  in  Spectral  Time- 

Domain  Simulations  of  Electromagnetic  Scattering.  Proc.  of  The  13’th  Annual  Review  of 
Progress  in  Applied  Computational  Electromagnetics,  Monterey,  CA.  pp.  926-933. 

•  D.  Gottlieb  and  P.  Fischer,  On  the  optimal  number  of  subdomains  for  hyperbolic  problems 

on  parallel  computers ,  Int.  J.  of  Supercomp.  Appl.  and  High  Perf.  Comp.  11(1997), 
pp.  65-76. 

•  D.  Gottlieb  and  S.  Abarbanel,  A  Mathematical  Analysis  of  PML  Methods.  J.  Comput.  Phys. 

134(1997),  pp. 357-363. 

•  M.  Carpenter,  J.  Norstrom,  adn  D.  Gottlieb,  A  Stable  and  Conservative  Interface  Treatment 

of  Arbitrary  Spatial  Accuracy ,  J.  Comput.  Phys.  148(1999),  pp. 341-365. 

•  S.  Abarbanel,  D.  Gottlieb  and  E.  Turkel,  Analysis  of  the  Error  for  Approximations  to 

Systems  of  Hyperbolic  Equations ,  J.  Comput.  Phys.  151(1999),  pp.  997-1007. 

111. 2  Submitted 

•  P.  G.  Dinesen  and  J.  S.  Fcsthaven,  2001,  A  Fast  and  Accurate  Boundary  Variation  Method 

for  Diffrative  Gratings  II.  The  Three-Dimensional  Vectorial  Case ,  J.  Opt.  Soc.  Am.  A  - 
submitted. 

•  J.  S.  Hesthaven  and  T.  Warburton,  2001,  High- Order /Spectral  Methods  on  Unstructured 

Grids.  I.  Time-Domain  Solution  of  Maxwell’s  Equations ,  J.  Comput.  Phys.  -  submitted. 

•  C.  H.  Teng,  A.  Ditkowski,  and  J.  S.  Hesthaven,  2000,  Modeling  Dielectric  Interfaces  in  the 

FDTD-Method:  A  Comparative  Study ,  IEEE.  Trans.  Antennas  Propaga.  -  submitted. 

•  P.  Dutta,  T.C.  Warburton,  and  A.  Beskok,  2000,  Analysis  of  Electroosmotically  Driven 

Micro-Channel  Flows ,  J.  of  Fluid  Mech.  -  submitted. 

•  A.  Beskok  and  T.C.  Warburton,  2000,  An  Unstructured  H/P  Finite  Element  Scheme  for 

Fluid  Flow  and  Heat  Transfer  in  Moving  Domains ,  J.  Comput.  Phys.  -  submitted. 

111. 3  Theses  Completed 

•  C.  H.  Teng,  Numerical  Methods  for  Wave  Problems  in  Complex  Geometries ,  PhD  thesis. 

Division  of  Applied  Mathematics,  Brown  University,  1998-2001. 

•  B.  Yang,  Spectral  Methods  and  Absorbing  Boundary  Conditions  for  Maxwell’s  Equations , 

PhD  thesis.  Division  of  Applied  Mathematics,  Brown  University,  1995-1998. 


15 


